COMPUTING ISOLATED SINGULAR SOLUTIONS OF 
POLYNOMIAL SYSTEMS: CASE OF BREADTH ONE* 
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Abstract. We present a symbolic- numeric method to refine an approximate isolated singular 
solution x = (xi of a polynomial system F = {/i , . . . , f n } when the Jacobian matrix of F 

evaluated at x has corank one approximately. Our new approach is based on the regularized Newton 
iteration and the computation of approximate Max Noether conditions satisfied at the approximate 
' singular solution. The size of matrices involved in our algorithm is bounded bynxn. The algorithm 

O^l ' converges quadratically if x is close to the isolated exact singular solution. 
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1. Introduction. 

Motivation and problem statement. Consider an ideal I generated by a 
■ polynomial system F — {/i , . . . , /„}, where f. k £ C[x±, . . . , x n ]. Suppose x = x e + x e , 

where x e denotes the isolated exact singular solution of F and x e denotes the error in 
the solution. The multiplicity p, of x e is defined as p — dim(C[x]/Q), where Q is the 
isolated primary component whose associate prime ideal is P = ( 
x n , e ), and the index p of x e is defined as the minimal nonnegative integer p such that 
P"'cQgSJ. ' 

In [36] [37], they compute the truncated coefficient matrix of the involutive system 
to the order p, and generate the multiplication matrices from its approximate null 
vectors. Then a basis of the approximate Max Noether space (Definition 12. ip of 
| / at x can be obtained from these vectors (Theorem 5.4 in [36]). Let y be the 

vector whose i-th element is the average of the trace of the multiplication matrix with 
respect to xi. In 1371, it has been proved that if the given approximation x satisfies 
qq . ||x — x e || = e, for a small positive number e, and the index p and the multiplicity 

\x are computed correctly, then the refined solution obtained by adding y to x will 
satisfy ||x + y — x e | = 0(e 2 ). Here and hereafter, || • || is denoted as the ^ 2 -norm. The 
size of these coefficient matrices in [37] is bounded by n( p ^™) x ( p ^ n ) which will be 
very big when p is large. Especially, when the corank of the Jacobian F' (x e ) is one, 
then p = p, which is also called the breadth one case in [4] [5]. 

In [18], we present a new algorithm which is based on Stetter's strategies [33] for 
computing a closed basis L = {L , . . . , of the approximate Max Noether space 

of / = (/i, ...,/„) at x incrementally in the breadth one case. The size of matrices we 
used in computing each order of Max Noether conditions is bounded by n x (n — 1), 
which does not depend on the multiplicity. Moreover, during the computation, we only 
need to store the input polynomial system F, the last n — 1 columns of the Jacobian 
F'(k) and the computed Max Noether conditions. Therefore, in the breadth one case, 
both storage space and execution time for computing a closed basis of the approximate 
Max Noether space are reduced significantly by the algorithm in [T5] . This motivates 
us to consider whether we can get rid of large coefficient matrices in |36l 137] and refine 
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approximate singular solutions more efficiently based on the computed Max Noether 
conditions. 

Main contribution. Suppose we are given an approximate singular solution x 
of a polynomial system F satisfying ||x — x e || = e, where the positive number e is 
small enough such that there are no other solutions of F nearby. We also assume that 
the corank of the Jacobian matrix F'(x e ) is one. In order to restore the quadratic 
convergence of the Newton method, we first apply one regularized Newton iteration (in 
Section [3. ip to obtain a new approximation x + y which also satisfies the assumptions 
above, and then compute the approximate null vector ri of the Jacobian F'(x + y) 
which gives a generalized Newton direction, and the step length 6 is obtained by 
solving a linear system formulated by the computed differential operators using the 
algorithm in [18] . We show that ||x + y + Sri — x e [| = 0(e 2 ). The size of matrices 
involved in our algorithm is bounded by n x n. The method has been implemented 
in Maple. Moreover, we also prove the conjecture in [5] that the breadth one depth- 
deflation always terminates at step fi — 1, where fi is the multiplicity. 

Structure of the paper. Section[2]is devoted to recall some notations and well- 
known facts. In Section [3] we describe an algorithm for refining approximate isolated 
singular solutions of polynomial systems in the breadth one case. Moreover, we prove 
that the algorithm converges quadratically if the approximate solution is close to the 
isolated exact singular solution. Some experiment results are given in Section [4] We 
mention some ongoing research in Section [5] 

2. Preliminaries. Let D(a) — D(ai, . . . ,a n ) ■ C[x] —5- C[x] denote the differ- 
ential operator defined by: 



D(ai,. . . , a n ) 



gaiH he 



ai l ■ ■ ■ a n l dx^ 1 ■■■dx^' 



for nonnegative integer array a = [a±, . . . ,a n ). We write 2) = {D(a), \a\ > 0} and 
denote by Span c (£)) the C- vector space generated by 2). Introducing a morphism on 
D that acts as "integral" : 

$-(D(a)):=l £> ( ai '---' Q! J _1 ' •••'"»)' ifa i> ' 
3 1 0, otherwise. 

As a counterpart of the anti-differentiation operator $ j , we define the differential 
operator <f j as 

*j(f(a)) := D(ai, . . . , aj + 1, . . . , a n ). 

Definition 2.1. Given a zero x e of an ideal I — (/i, . . . , f n ), we define the Max 
Noether space [23] of I at x e as 

A 4e (/) := {L G Span c (S))|i(/) x= x c = 0, V/ £ /}. 

Conditions equivalent to L(f) x —^ e = 0, VL G A% e (I) are called Max Noether condi- 
tions [23 . The space {L^ c \ L G Ai t (/)}, where Ls e (/) := £(/) x= a e , is also called 
the dual space of I at x e [HI [20] [22] EH [33] • For a nonnegative integer k, Al^(I) 
is a subspace of A Xe (/) which consists of differential operators with differential order 
bounded by k. Obviously, A^ 0) (7) = D(0, . . . , 0). We have that 



dim c (A Xe (/)) = fi, 



(2.1) 
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where [i is the multiplicity of the zero x e . 

Definition 2.2. J7H ^4 subs-pace A o/Span c (£>) is said to be closed if and 
only if its dimension is finite and 

L G A => $,•(£) G A, j = l,...,n. 

Suppose Span(L 5 Lx, • ■ ■ , Lu—x) IS closed and L , . . . , L^-x are linearly independent 
differential operators which satisfy that £i(/j) x =x e = 0, j = 1, . . . , n, i = 0, . . . , (J, — 1, 
then due to the closedness, Lj(g • /j) x =* e =0, Vq £ C[xi, . . . , £„]. Hence, A ie (I) = 
Span(L ,ii, . . . 

Lemma 2.3. Let F'(x. e ) be the Jacobian of a polynomial system F = {fx, . . . , /„} 
evaluated at x e . Suppose the corank of F' (x e ) is one, i.e., </ie dimension of its 
null space is one, then dim(Al (/)) = dim(A^ fc (I)) + 1 f or I < k < fi — I and 

dim(A^ (7)) = dim(A^ 1 1 (/)), for k > \x. Hence we have \i = p. 

Proof. Lemma 12.31 is an immediate consequence of [32] Theorem 2.2] and [SJ 
Lemma 1]. □ 

Theorem 2.4. Jjjjjj Suppose we are given an isolated multiple root x e of the 
polynomial system F — {/l, ••-,/«} with the multiplicity fi and the corank of the 
Jacobian F' (x e ) is one, and Lx = D(1,0,...,0) G A^J). We can construct the 
k-th order Max Noether condition retaining the closedness incrementally for k from 2 
to [1 — 1 by the following formulas: 

L k =P k + a k , 2 D(0, 1, . . . , 0) + • • • + a k , n D{0, . . . , 1), (2.2) 

where P k has no free parameters and is obtained from previous computed Lx, . . . , Lk-i 
by the following formula: 

(2.3) 

where 

$i(P fe ) = L fc _!, Q fc j = $j(P fc ) = a 2 ,jL k - 2 + ■■■ + a k -x,jLx, 2 < j < n. (2.4) 

Here ^j(Qk.j)a 1 =---=a j _ 1 =o means that we only pick up differential operators D(a) 
in Qk,j where ax = • • ■ = cx.j—x = 0. The parameters a k j, j = 2, . . . ,n are deter- 
mined by checking whether [-Pfe(/i) x =x e) ■ • ■ , Pk(fn)x=Z e ] T can ^e written as a linear 
combination of the last n — 1 linearly independent columns of F'(x e ). 

Suppose x is an approximation of x e and ||x — x e || = e <C 1, we can use the 
algorithm MultiplicityStructureBreadthOneNumeric in [18 to compute a closed basis 
{Lq, . . . , Ln-x} of the approximate Max Noether space of I at x. Since the errors in 
the matrix of the linear system 



dF(i) 5F(x) 

"h\* Jx= 



[1, a k ,2, ■ ■ ■ ,ak, n ] T = 0, 



dx 2 ' dx n 

used in Theorem 12.41 are bounded by O(e) and 

L k = P k + a k , 2 D{0, 1, 0, . . . , 0) + • • • + a k , n D(0, ...,0,1), 

is determined by its right singular vector [1, a k ,x, ■ ■ ■ , Ok,n] T corresponding to its small- 
est singular value, we have 

\\Lk(F)x=*\\ = O(e), (2.5) 
according to [TTJ Corollary 8.6.2]. 
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3. An Algorithm for Refining Approximate Singular Solutions. Suppose 
we are given an approximate solution 



where x e denotes the error in the solution and x e denotes the exact solution of the 
polynomial system F — ...,/„} with the multiplicity p and the index p. In this 
section, we present a new method to refine x in the breadth one case, i.e., p = p. 

Let A — F'(Sc) be the Jacobian matrix of F evaluated at x and b = — F(x). 
Suppose the error in the solution is small enough, i.e., ||x — x e || = e -C 1, and A is 
invertible, then Newton's iteration computes 



and ||x + y — x e || = 0(e 2 ) according to the well-known Kantorovich theorem [13] . 
However, if A is singular, then the convergence of Newton iterations is linear rather 
than quadratic. 

Rail [27] studied the convergence properties of Newton's method at singular 
points. Some modifications of Newton's method to restore quadratic convergence 
have also been proposed in [U El El El H2J H31 HH EH EH EH 1^1 EI] . In [13], a bor- 
dered system was introduced to restore the quadratic convergence of Newton's method 
when A has corank one approximately and x is a simple singular solution. It is clear 
to see that the regularity condition in [13] can not be satisfied if the multiplicity is 
larger than 2. 

For simplicity, we make an assumption throughout this section. 

Assumption 1. Suppose we are given an approximate singular solution x of a 
polynomial system F satisfying ||x — x e || = e, where the positive number e is small 
enough such that there are no other solutions of F nearby. Moreover, we assume that 
the corank of the Jacobian matrix F'(x e ) is one. 

Let A = F'(x) be the Jacobian matrix of F evaluated at x and its singular values 
be CTi, . . . , a s . Under Assumption[TJ we have ||F(x)|| = 0(e), at = 0(1), 1 < i < n—1 
and a n = 0(e). 

Remark 3.1. Notice that the notation 0(g) denotes that the value is bounded 
above by g up to a constant factor, while 0(g) denotes that the value is bounded both 
above and below by g up to constant factors. 

3.1. Regularized Newton iteration. Under Assumption [l] F'(x) is approxi- 
mately singular. Instead of using (|3.ip to compute y, we apply Tikhonov regulariza- 
tion |34j to solve the minimization problem 



to obtain y, where A = F'(ic) and b = — F(x). The real number A > is called the 
regularization parameter. 

Theorem 3.2 (Regularized Newton Iteration). Under Assumption^ if we choose 
the smallest singular value a n of F'(it) as the regularization parameter, the solution 
y of the following regularized least squares problem 



x = x e + x, 



y = A- 1 b, 



(3.1) 



min Py-b|| 2 + A||y|| 2 



(A* A + a n I n )y = A*b 



(3.2) 



satisfies 



||y||=0(e), ||F(x + y)||=0(e 2 ), 



(3.3) 
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where A* is the Hermitian (conjugate) transpose of A = F'(x), I n is the nxn identity 
matrix and b = — F(x). 

Proof. Suppose A = U ■ £ • V* is the singular value decomposition of A where 
£ = diag{cri, . . . , cr„}, then the solution of (|3 . 2() is 

y = V-(Y?+a n I n )- 1 -^-U* • b. (3.4) 

Since er^ = 6(1), l<i<n — 1, a n — 0(e) and ||b|| = 0(e), we have 



where b = . . . , b n ] T = U*b and ||b|| = ||b|| = 0(e). Hence, ||y|| = 0(e). 
From the Taylor expansion of F at x, we have 

J F(x e ) = -b + A(x e -x) + 0(e 2 ). 

Hence 

||-b + A(x e -x)|| = 0(e 2 ). 

Furthermore, we have 

||-Erb + £.V*(x e -x)|| =0(e 2 ). 

Since a n — 0(e) and ||F*(x — x e )|| = ||x — x e || = e, we derive that the last component 
of the vector b satisfies 



IM = 0(e 2 ). (3.5) 



Since 



we have 



Ay - b = U • diag ^ -g- , . . . , -j- } ■ b, 



||Ay-b|| 2 = ]T 



n / \ 2 

r-» / o- n \bi\ \ 



where 



07 + a>, 



0(e), for i = 1, . . . ,n — 1, 



and ||b|| = 0(e). Although 



= 9(1), 



we have from (|3.5p that |6 n | = 0(e 2 ), hence 

||Ay-b||=0(e 2 ). (3.6) 



G 



Finally, from the Taylor expansion of F at x, we have 

||F(x + y)|| < ||-b + Ay||+0(e 2 ) = 0( £ 2 ). 

□ 

According to Theorem 13.21 after applying one regularized Newton iteration to 
F and x, we get y satisfies (|3.3|) . and the new approximate singular solution x + y 
satisfies 



If 



|x + y-x e || < ||x-x e || + ||y|| =e + 0(e). 



|x + y-x e || = 0(s 2 ), 



then we have already achieved the quadratic convergence. However, the convergence 
rate of the regularized Newton iteration is linear too when the Jacobian matrix is 
near singular. Hence, in most cases, we will have 



|x + y-xe||=0(e). 



(3.7) 



We show below how to restore the quadratic convergence when the computed approx- 
imate singular solution x + y satisfies (|3.3[) and (|3.7p . 

If L\ 6 is not £>(1,0, ... ,0), as pointed out by Stetter in [33], we can 

compute the right singular vector of F' (x + y) corresponding to its smallest singular 
value a' n , denoted by ri satisfying ||ri|| = 1 and 

HF'Cx + yJnll =a' n = 0(e). (3.8) 

Let us form a unitary matrix R — [r±, . . . , r n ] and perform the linear transforma- 



tion 



It is clear that 



H(z) = F(Rz). 



z e = R 1 x e 



is an exact root of H(z) and 

z = i?- 1 (x + y) 

is an approximate root of H(z). Moreover, we have 

||z-z e || = Hfl-^x + y-xe)!! = ||x + y-x e || =0(e), 



(3.9) 



(3.10) 



(3.11) 



(3.12) 



|tf(z)|| = ||F(x + y)|HO( £ 2 ), 



and 



dH{z) 



\F\k + y)v l \\=a' n = 0{e). 



(3.13) 



(3.14) 
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Hence, the condition (|3.7[) is equivalent to (|3.12l) . Here and hereafter, we always 
assume that z satisfies 

||z-z e || = 6(e). (3.15) 

Theorem 3.3. The root z e defined in A3.10\) is an isolated singular solution of 
H with the multiplicity /i and the corank of H'(z e ) is one. 

Proof. Since H'{z e ) = F' (x e )i? and R is a unitary matrix, we derive that the 
corank of H'(z e ) is one. Let fj! be the multiplicity of z e , and {L , L\, . . . , Lu'-i} be 
a closed basis of the Max Noether space of H at z e . The operator T R : Span c (£>) — » 
Span c (2)) is defined by: 



T R (D(a)) := L fl 



Q<*l-\ fOn 



ai! • • • a„! dz" 1 ■ ■ ■ dz„ 
ai! • ■ • a n \ d(r{ • x) Ql • • • <9(r* • x) Q " 

1 ^ ^l+-+ffn 



ai! • • • a 



I0I=M 



where c« is the coefficient of — in the expansion of », , 3 1+ * " — _ - . Since 
i?(z) = F(Rz) and x — Rz, according to multivariate chain rules, we have 

r fl (i fe )(^)x=x e = L k (H) z=ic = 0, 

and for 1 < i < n, 

%(r fl (L fc )) =r fl ^r w *i(L fc )J 

/fc-2 \ 

=Fr I 2j(afe-»,2r2,j H h ak-i,nr n ,j)Li +ri tj Lk-i J 

fc-i 

= 2j(afc-i,2'*2,j H 1- ak-i,nr n ,j)TR(Li) + rijT R (L k _i), 



where < < //' — 1. Hence, {T R (L ), T r (Li), . . . , r^L^z-i)} is a closed basis of 
(L) and // < /i. On the other hand, since L(x) = H(R^ 1 x), we derive that 
fi < //. Hence, // = /i. □ 

Remark 3.4. Since H'{z) = F' (x+y)i? and R is a unitary matrix, we derive that 
the singular values of H'(z) are the same as those o/-F'(x+y) and the corank of H'(z) 
is one approximately. Suppose {Lo, Li, ... , L^_i} is a closed basis of the approximate 
Max Noether space of H at z, where Lq = L>(0, . . . , 0) and L\ = D(l, 0, . . . , 0). From 
the proof of Theorem \3.3\ and &2.5\) , we have 

Lfl(L fe )(F) x=jc+ y = L k (H) z=z = 0(e), 

and 

%(r R (L fe )) G Span{r R (L 1 ), . . . , T R (L k ^)}, 
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where < k < fi — 1 and 1 < j < n. Hence, {Tji(L ), T^(Li), . . . ,Tji(L^i)} is a 
closed basis of A x+y (I) . 

Remark 3.5. It should be noticed that Theorem \S.3\ holds as long as R is a 
regular matrix. However, if we choose a unitary matrix R, then it is much easier to 
compute the inverse of R since R^ 1 = R* . 

It is interesting to notice that, after running one regularized Newton iteration, 
the last n — 1 elements of the solution z have already been refined quadratically. 

Theorem 3.6. Suppose z e and z are defined in \3.1U\) and \3.11\) respectively. 
Under Assumption^ we have 



pi, 



0(e), 



(3.16) 



\k,e - z%\ = 0(e 2 ), for i = 2, . . . ,n. 
Proof. From the Taylor expansion of H(z) at z, we have 
H(z e ) = H{z) + H\z)(z e - z) + 0(e 2 ). 
Since H{z e ) = and \\H(z)\\ = 0(e 2 ), we have 

||i/'(z)(z e -z)||=0( £ 2 ). 
From (I3.14|) and (|3.15p . we have 

dH{z) 



and 



dzi 



dH(z) dH(z) 



{Zl,e - h) 



8Z2 



dz n 



[Z2,< 



According to Remark 13.41 the matrix 
that (13.17)) is correct. The equation ^£ 



Z2, 



8H(z) 

dz2 



(3.17) 



3H(z) 



0(£ 2 ). 



is of full column rank, so 



U) follows from (j3"T5l) and (|3~T7| . □ 



If the multiplicity /i is larger than 2, the regularity assumption in p~3] will not be 
satisfied. The violation of the regularity assumption is caused by the existence of the 
higher order Max Noether condition. It is interesting to notice that the left singular 
vector of the Jacobian matrix H' '(z) corresponding to the smallest singular value can 
be used to prove the following theorem. 

Theorem 3.7. // the multiplicity of the singular root is larger than 2, under 
Assumption^ we have 



\Li(H) 



dH{z) 

dzi 



0( £ 2 ). 



Proof. If fi > 2, according to Theorem 
Max Noether condition such that 



and (12.51) . there exists a second order 



\L 2 (H) z=i 



2dl 



a-2,2 



d 
dz 2 



a 2: „^- ) (H) z=i 



= 0(s) 



(3.18) 



Suppose u n is the left singular vector of H'(z) corresponding to the smallest 
singular value a' n and ||u n || = 1, then 



dH(z) 



dzi 



0(e), 1 < i < n. 



From (|3~T8| and (f37T9|) . we have 



3 2 H{i) 



dzf 



0(e). 



Therefore, we get 



d 2 H{i) 

■ dz\ 



(zi,e - ZiY 



0(e 3 ). 



From the Taylor expansion of H (z) at z, we have 

H(z e ) = H(z) + H'(z)(z e - z) + H"{z)(z e - z) 2 + 0(e 3 ), 



(3.19) 



(3.20) 



(3.21) 



(3.22) 



where (z e — z) 2 denotes the vector of all monomials with degree 2 and H"(i) consists 
of all second order derivatives of H evaluated at z. 

According to Theorem l3.61 all elements in (z e — z) 2 are 0(e 3 ) except the first one. 
Combining with p.21j) . we have 

|<i/"(z)(z e -z) 2 | =0(e 3 ). (3.23) 

On the other hand, the Taylor expansion of H(z) at z e shows that 

H(z) = H(z e ) + H'(z e )(z - z e ) + H"(z e )(z - z e f + 0(e 3 ). 

Since the corank of H'(z e ) is one, suppose u e is the left null vector of H'(z e ) and 
llu.ll = 1, then 



u* e H'(z e ) = 0. 



Notice that 



Kif'( Z )|| = \\u* e [H'(z) - H'(z e )}\\ < \\H\z) - H'(z e )\\ = 0(e), 

and H'{z) has corank one approximately, so that ||u„ — u e j = 0(e). Moreover, using 
the same analysis above, we obtain that 

|u:tf"(z e )(z-z e ) 2 | = 0(e 3 ). 

Hence, we have |u*If(z)| = 0(e 3 ). Noticing ||if(z)|| = 0(e 2 ), we get 

\KH(z)\ < |(u„ - u e )*H(z)\ + \u* e H(z)\ = 0(e 3 ). 

Combining |3~22l) . (|3~23)) and (j3~24| . we have 

|<tf'(z)(z e -z)| = 0(e 3 ), 

which is equivalent to 

K<(z e -z)| = 0(e 3 ), 



(3.24) 
(3.25) 
(3.26) 
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where v„ = [1,0, ...,0] T is the right singular vector of H'(z) corresponding to a' n . 
Hence, \a' n (z 1 „ e — Z\)\ = 0(e 3 ). Based on (|3.16[) . we have 

< = 0(£ 2 ). (3.27) 



Moreover, from (|3.14[) . we have 

= 0(e 2 ). (3.28) 
□ 

It is amazing to notice that not only the first order Max Noether condition com- 
puted according to Theorem l2. 41 satisfies (|3.28|) . but also all other Max Noether con- 
ditions up to the order fi — 2 > satisfy similar conditions: 

\\U(H) z=i \\ = 0(e 2 ), fori = 0,...,/*-2. (3.29) 



\Li(H) z=± \\ = 



dH(z) 



3.2. An Augmented Polynomial System. To prove (|3.29p inductively, we 
need to introduce an augmented polynomial system and prove the following theorem. 

Theorem 3.8. Let us assume that H(z) is a polynomial system which has z e 
as an isolated exact singular solution with the multiplicity /i, the corank of H'(z e ) 
is one. Let L be the ideal generated by polynomials in H and {Lq, Li, . . . , i} 
be a closed basis of Aj E (7), where Lq = D(0, . . . , 0), L% = D(l, 0, . . . , 0) and Lk = 
Pk + flfc.2-D(0, 1, . . . , 0) + • • • + afc j71 D(0, . . . , 1) constructed according to Theorem \2.4\ 
The augmented polynomial system 

H(z), 

G(z,A) := { H'(z) ■ A, (3.30) 
Ai - 1, 

where A = [Ai, . . . , A„] T has an isolated singular solution (z e , A e ) with the multiplicity 
/i — 1, where A e = [1, 0, . . . , 0] T . // fi > 3 then the Jacobian matrix G'(z e ,A e ) has 
corank one and 

d d d , 

Li = tt- +2a 2 ,2 7n- + --- + 2a 2,n7TT- (3.31) 
ozi d\ 2 o\ n 

satisfies Li(G)^ z \)=(z a ) = 0- Moreover, starting from Lq = D(0, . . . , 0) and L\, for 
2 < k < /i — 2, the k-th order Max Noether condition of G at (z e , A e ) retaining the 
closedness has the following form: 

d d d d 

Lk = Pk + ak,2s 1 \~ak,nT{ Kfc + l)a fc+ i i2 7rT-H H(A;+l)afc + i, n — — (3.32) 

oz 2 oz n d\ 2 o\ n 

where 

Pk = Pk + ^n+2(Qk,n+2) + * n+z{Q fc,n+3)a„ + 9=0 + ' " ' + ^2n{Q k,2n)a n + 2 = - ■ -=a 2n -i=0 

(3.33) 

and 



Qk,n+j = $n+j(Pk) = 2a 2 ,jL k _i H h ka k jLi, 2 < j < n. 



(3.34) 
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Proof. The Jacobian matrix of G(z, A) at (z e , A e ) is 
G'(z e , A e ) 

where H"{z e )-\ e = 



H'%) ' 

H"(m e ) ■ A e H'(z e ) 

XT 



a-=g(z e ) 

5 c*2:i<9.s 71 



D(l, 0, . . . , 0) £ A^ 1 ' (J), the first column of H'(z e ) is a zero vector and the remaining 
columns of i?'(z e ) are linearly independent. Moreover, since \l = [1,0,... ,0], the 
last In — 1 columns of G'(z e , A e ) are linearly independent and its corank is less than 
one. 

If /i > 3, the second order Max Noether condition of H at z e has the form L 2 = 
D(2, 0, . . . , 0) + 02,2^(0, 1, 0, . . . , 0) + • • • + a 2 ,„-D(0, ...,0,1). From L 2 {H) z=ie = 0, 
we have 



Since the corank of H'(z e ) is one and L\ 



l d 2 H(z e ) dH(z e ) 
2 dz\ 2,2 dz 2 



, dH(z e ) 
+ a 2 . n — q = 0. 



The vector v = [1, 0, . . . , 0, 202,2, ■ • ■ , 2a 2yn ] T is a null vector of G'(z e , A e ). Therefore, 
the Jacobian matrix G'(z e , A e ) has corank one and the first order differential operator 
L\ in (|3.3ip satisfies 



Li(ff'(z)-A) 



(z,A) = (z c ,A e ) 



2L 2 (H) z= z e = 0. 



Hence, we have 



Li{G) 



(z,A) = (z e ,A e ) 



0. 



(3.35) 



(3.36) 



Using similar arguments in [18 for proving Theorem l2.4l we can show that the dif- 
ferential operators L k defined by formulas (|3. 321) . p.33p and (|3.34p retain the closed- 
ness. It should also be noticed that L k always contains the differential monomial 
D(k, 0, ...,0) and there are no differential monomials D(i,0, ...,0) for i < k con- 
tained in L k . Otherwise, we can reduce them by Li. Moreover, is not contained 

in any L k , otherwise, L k (\x), z x)=(z a ) ^ 0- Hence, due to the closedness, there 
are no differential operators D{a\, . . . , a n , a„+i, . . . , cx2 n ) with a n+1 > contained in 
any L k . 

Now let us show that the constructed differential operators L k satisfy 



L k (G) 



(z,A) = (z«,A«) 



= 0, for 1 < k < fi - 2. 



(3.37) 



From (|3.36[) . we can see that (|3.37p is true for k = 1. Moreover, it is easy to check 
that 



L k (H) 



(z,A) = (z e ,A.) 



= Pfe + a k - 



_8_ 

dz 2 



ak,-. 



d 



(#%=z e = 0, 



(3.38) 



and 



ife(Ai) 



(z,A) = (z e ,A e ) 



(3.39) 
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Based on formulas (|3.32[) . (|3.33[) and (|3.34l) . we have 
L k (H'(z) ■ A) (ZiA)=( - e Xe) = 

L fc — +^(2a 2 ,ji&-i H h fca fej ii + (fc + 1) afe+ij)— J (iJ) 2 



Let us set 



d d 
Q fc+ i =L k — + y\{2a 2 ,jL k _ 1 + --- + ka kd L 1 )— . (3.40) 

j=2 3 



We show (Proposition 15. II in Appendix) that 

Q k+1 = (k + l)P k+1 . (3.41) 

Hence, we have 

L k (H'(z).\)^ x)={ , eX) 

= ((fc + l)P k+1 + (k + l)a W7 ^ + -.. + (k + ljofc+i,,^) (fl)-=*. ( 3 - 42 ) 
= (fc + l)i fe+ i(JJ) z=2e =0. 

From (|3~5gj) . (j3~31?)) and (|3~l2"j) . we derive that (1337)) is true forl<fc</x-2.0 

Corollary 3.9. Suppose F(x) is a polynomial system which has x e as an iso- 
lated exact singular solution with the multiplicity fi and the corank of F' (x e ) is one. 
Let Y\ be the null vector of F' '(x e ) and \\v\\\ = 1. For any random vector h G C" 
satisfying h*ri ^ 0, the augmented polynomial system 

r f(x), 

J(x,i/) := { (x)-i/, (3.43) 

/las (x e , h ^ ) as an isolated singular solution with the multiplicity jJL — 1. 

Proof. Let {r 1; . . . , r„} be a normal orthogonal basis of C™, then h = (h*rx)ri + 
• • • + (h*r n )r n . If h*ri ^ 0, performing the linear transformation 



ri „ _ h*r 2 r r _ h*r„ 

h*n' 12 h'n 1 h*ri 1 



where R = 
augmented polynomial system 



x = Rz, v = RX, 

is a regular matrix, we obtain the 



H(z), 

G(z,A) := { H'(z) ■ A, 
. Ai-1, 

where 

H(z)=F(Rz), H'(z) ■ A = F'(x) • R ■ R~ 1 v 1 \ x - 1 = h*u - 1. 

According to Theorem 13. 8[ we know that (z e ,A e ) is an isolated singular solution of 
G with the multiplicity \i — 1, where z e = .R _1 x e and A e = [1, 0, . . . , 0] T . Hence, by 
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Theorem 13.31 and Remark 13.51 ^x e , jpr^) is an isolated singular solution of J{x.,v) 
with the multiplicity fi — 1. □ 

Remark 3.10. It is well known that the augmented polynomial system J(x,u) 

defined in \3.43\j has an isolated singular solution ^x e , pr^-^ with the multiplicity less 

than \i, see \VI\ [3]/. Here, we proved the conjecture in '51 that the multiplicity of the 
singular solution of the augmented polynomial system \3.43ty drops by one exactly in 
the breadth one case. 

Remark 3.11. For the system H(z) and its approximate singular solution z de- 
fined in &3. 9\) and &3.11\) , the augmented polynomial system defined in h3. 30\) has (z, A) 
(X = [1, 0, . . . , 0} T ) as an approximate solution. Suppose {Lo, . . . , L^-\} is a closed ba- 
sis of the approximate Max Noether space of the system H atz constructed according to 
Theorem \2.4\ from Lq = D(0, . . . , 0) and L\ = D(l, 0, . . . , 0), then {Lo, Li, ... , £^-2} 
constructed according to Theorem \3.8\ is a closed basis of the approximate Max Noether 
space of the system G at (z, A), satisfying 

\\Lk{H) { ^ x)=(i:x) \\ = \\L k {H) z=i \\ = 0{el 
ll^(tf'%,A)=(*,A)ll = \\(k + l)L k+1 (H) z= 4 = 0(e), 
l|£*(Ai)(. ( A)=(*,A)ll=0, 

for 1 < k < (j, - 2. 

Theorem 3.12. Let F(x) be a polynomial system which has x e as an isolated 
exact singular solution with the multiplicity \x and the breadth one. Suppose x is an 
approximate solution of F which satisfies 

[[x - x e [[ = 9(e) and ||F(x)[| = 0(e 2 ), (3.44) 

for a small positive number e. Let <j\, . . . ,a n be the singular values of F'(x) satisfying 
o~i = 6(1), I < i < n — 1 and a n — O(e). Suppose i"i is the right singular vector 
corresponding to a n . We form a unitary matrix R = [ri,...,r„] and set H(z) = 
F(Rz). Suppose {Lq, . . . , L^i} is a closed basis of the approximate Max Noether 
space of the system LI at z = i? _1 x constructed according to Theorem \2.4\ from 
L = D(0, . . . , 0) and L x = D(l, 0, . . . , 0), then 

||i<(-E0-=s|| = °( £2 )> for i — 0, . . . , fi 2. 

Remark 3.13. Under Assumption^ according to Theorem \3.2l we can always 
perform the regularized Newton iteration to obtain an approximate singular solution 
x satisfying \3. 44\ )- Moreover, it should also be noticed that all discussions in Section 
\3.1\ after Theorem \3.2\ are valid if we start with an approximate singular solution 
satisfying \3.44^ - 

Proof. According to (|3.13[) and Theorem 13.71 we know that Theorem 13.121 is true 
for fi = 2 and /i = 3. 

Now let us assume that Theorem 13. 121 is true for fi = k and k > 3. For fj, = k + 1, 
we form the augmented polynomial system G{z, A) = {H(z), H'(z) • A, Ai — 1}. 

According to Theorem 13.31 the root z e defined in (|3. 10[) is an exact singular 
solution of H(z) with the multiplicity /j, and the corank of H'(z e ) is one. Let v be 
the null vector of H'(z e ) and ||v|| = 1. Since 

||ff'(i)v|| = \\[H'(z) H'(z e )]y\\ = O(e), (3.45) 
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dH{i) 
dz 2 : 



dH{i) 



is of full column rank, combining with p,14[) . we derive that 
vi = 6(1), and v, = 0(e), for 2 < i < n. (3.46) 



Set h= [1,0,..., 0] T , we have h*v = i> x =6(1) ^ 0. 

According to Corollarv l3.91 the augmented polynomial system G(z, A) has (z e , 
(z e , A e ), where 



X, = 



i "Hi "Hhl 

1, , . . . , 

Vl Vi 



as an isolated singular solution with the multiplicity fi — 1, which is equal to k. 
According to Remark 13.111 (z,A) (A = [1, 0, . . . , 0] T ) is an approximate solution of 
G(z, A). Moreover, by (|3~T5|) and (pHrJ|) . we have 



||(z,A)-(z e ,A c )|h 
Furthermore, from (|3.13p and Q3.28[) . we have 



+ IIA - A e 



6(e). 



||G(z,A)|| = W||tf(z)| 



dH(z) 



0(e 2 



We have assumed that Theorem 13.121 is true when the multiplicity is equal to 
k. Therefore, for the augmented polynomial system G(z,A), we can form a unitary 
matrix R with = -[1, 0, . . . , 0, 2a2,2i ■ • ■ > 2a2 ;n ] T as its first column, where a 



1 + 4(a| 2 + • • • + 0% n ), then generating a new system J(w) = G(R w) which has an 
approximate singular solution w with the multiplicity k. By the inductive assumption, 
we have 

||Li(J)w=w|| = 0(e 2 ), for < i < k - 2, 

where Li is the i-th. Max Noether condition of J at w constructed by Theorem 12.41 
from Z = D(0, . . . , 0) and L x = D(l, 0, . . . , 0). According to TheoremE31 



(«,A)=(fi,A)- 



Since {Lo, ^l, ■ • ■ , Lk-i} and {r^(Lo), T^(Li), . . . , r^(Lfc-i)} are both closed basis 
of the approximate Max Noether space of the system G at (z, A), and 



1 - 



Fx(L ) = L , and T R (Li) = -Li, 

we derive that r^(L^) is a linear combination of {Lg, L\, . . . , Li} (Proposition 15.21 in 
Appendix). Hence, we have ||£,(G), „y)=(z aJ = ^(e 2 ), and 



(z,A) = (z,A) 
1 



Li(H'X) 



(z,A) = (z,A) 



0(, 2 ). 



Therefore, Theorem 13. 121 is true for fi = k + 1. □ 
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3.3. An Algorithm for Refining Approximate Singular Solutions. 

Algorithm 1. MultipleRootRefinerBreadthOne 
Input: An approximate solution x of a polynomial system F which is close to an 
isolated exact singular solution of F with the multiplicity (i in the breadth one case, 
and a tolerance T. 
Output: Refined solution x. 

1. Regularized Newton Iteration: Solve the regularized least sguares problem 

(A*A + a n I n )y = A*b, 

where b = — F(x), A* is the Hermitian (conjugate) transpose of A = F'(i), 
I n is the n x n identity matrix and o~ n is the smallest singular value of A. 

2. Compute the null vector T\ of F' (x + y) with respect to t, form a unitary 
matrix R with V\ as its first column and perform the linear transformation 

H(z) := F{Rz), 

and set z :— i? _1 (x + y). 

3. Construct a closed basis of the approximate Max Noether space of I = (hi, . . . , h n) 
at z with respect to t: 

A^" X) (/) ^Span^cii,---,^-!) 

by Algorithm MultiplicityStructureBreadthOneNumeric in \18^ . 
4- Solve the linear system 



v = -L M _i(fO-=i. (3-47) 



dz 2 ' " dz n 

where v = [v\, ■ ■ ■ , v n ] T and is the differential operator of order \i com- 
puted by formulas in Theorem \2.4\ Set 6 := — . 
5. Return 

x := x + y + 5t\. 

Remark 3.14. The size of matrices involved in the algorithm MultipleRootRe- 
finerBreadthOne is bounded by nxn, whereas the size of matrices used in the deflation 
method is bounded by (u n) x (// n) |T7| /. 

Remark 3.15. In fact, in order to keep the sparse structure of the original poly- 
nomial system, we should avoid performing the linear transformation. Moreover, it is 
expensive to compute and store all Max Noether conditions. Since we only need their 
evaluations to solve Jff.^7| ) ; it's possible to compute and store only the necessary eval- 
uations of these Max Noether conditions. We will discuss these issues in forthcoming 
papers. 

3.4. Quadratic Convergence of the Algorithm. 

Theorem 3.16. Under Assumptions^ the refined singular solution x returned 
by Algorithm MultipleRootRefinerBreadthOne satisfies 



||x-x e ||=0( £ 2 ). 



(3.48) 
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Proof. According to Theorem 13. 121 we have Li(H) z=z = 0(e 2 ), for < i < \i — 2. 
Since 

Ofc(Li) G Span(L , ■ • ■ , ^-2), for 1 < fc < n, 

we have 

||Li((«*-JSk)H).=i|| - \\M L i)( H )*=4 = 0(e 2 ), for 1 < k < n, < * < (3.49) 

The matrix in (|3.47[) is of full rank. We solve the linear system (|3.47l) to obtain 
the vector v = [vi, . . . , v n ] T such that £ M (-ff) z=z = for 

d 

:= L^-i + vi ■ + v 2 ■ 7^- -\ h 

It should be noticed that the vector v satisfies ||v|| = 0(e) since \\L fJl -i(H) z=z \ 
0(e). Moreover, 



(3.50) 



$fc(L M ) G Span(L , ■ • ■ , i/t-a), for 2 < fc < n, 



we have 



\\Lpdzk - z k )H), =i \\ = \\<t>k(L^)(H) z=z \\ = 0(e 2 ), for 2 < k < n. 
For k = 1, since ||$i(wiP M )( J ff) z =z|| = ||wi£ M -i(-ff)z=z|| = 0(e 2 ), we have 

II =0(e 2 ). 



(3.51) 
(3.52) 



\\L,(( Zl - Zl )H) z=z \\ = W^iL^H)^ 
From (pH9"l) and (j3~52"1) . for i = 0, 1, . . . , /i - 2, ,u, we have 

WHp ■ H) z=z \\ = 0(e 2 ), Vp G {(zj - z x )^ ■ • • (z„ - z n ) Q ", ax > 0, . . . , a„ > 0}. 
Especially, we have 

\\M ll+1 -L l {v{z),) z=z \\=0{E 2 ), 

where A/ M +i is the coefficient matrix of the Taylor expansion of the system H and all 
its prolongations up to the degree /i at z, and 

v(z) /1 = [(zi - ii) M , (zi - zi) Al_1 (z 2 - z 2 ), • ■ • ,zi - 2i, . . . ,z n - z n , l] T . 

It is important to notice that, based on the closedness conditions, we obtain the 
null space of M^ + \ with matrices of size n x n instead of generating the big matrix 
M^ + i. Similarly to the analysis in [37, Remark 18], the trace of the multiplication 
matrix M Zl formed from approximate null vectors -Li(v(z) At ) z=z has the following 
property 



-Tr(M zl ) = -Tr(M 21 ) + 0(e 2 ) = -z x ,< + 0(e 2 ). 



(3.53) 



It is interesting to notice that, by using the approximate basis {Lq, . . . , L^} and 

the normal set |l, . . . , J^-i 1> we can form the multiplication matrix 

Oh ■•• 



M z 



lo 






h 



















h 
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where U is the coefficient of in Li. Hence, the trace of M Zl is v\. Therefore, there 
is no need to form the multiplication matrix! According to (|3 . 53|) . we have 

-+i l£ = 0(e 2 ). (3.54) 
M 

Since the last n — 1 elements of z have already been refined quadratically, by 
updating Z\ := Z\ + 5 for 5 :— we have 

||x-x e || = \\R(z-z e )\\ = 0(e 2 ). 

□ 

Remark 3.17. TTie algorithm MultipleRootRefinerBreadthOne also works well for 
some over determined 'polynomial systems, i.e., the number of polynomials is bigger 
than the number of variables, see the last example Menzell in Table [^J 

4. Examples. The following experiments are done in Maple 13 under Linux for 
Digits :— 15. Let t and s be the number of polynomials and variables respectively, 
H be the multiplicity. Systems DZ3, Dayton2 and DLZ are quoted from [3J HI [5], 
Menzell and SY5 are cited from [3T] and [3T] respectively. Other examples are cited 
from the PHCpack demos by Jan Verschelde. 



System 


Zero 


t 


s 


M 


# Digits 


Ojikal 


(1,2) 


2 


2 


3 


2 -> 5 -> 11 -> 15 


Ojika2 


(1,0,0) 


3 


3 


2 


2 -> 5 -4- 10 -4- 14 


Ojika3 


(-2.5,2.5,1) 


3 


3 


2 


2 -> 4 -> 9 ->■ 14 


Ojika4 


(0,0,10) 


3 


3 


3 


2 -> 3 -> 7 -> 13 


Decker2 


(0,0) 


2 


2 


4 


2^5^15 


DLZ 


(0,0) 


2 


2 


10 


2 ->■ 5 -> 16 


DZ3 


(2Vl , VS \/7 , 2V5\ 


2 


2 


5 


2 5 13 


Dayton2 


(0,0,0) 


3 


3 


5 


2 3 -> 7 13 


SY5 


(1.1) 


2 


2 


2 


2 -> 5 ->■ 11 -s- 14 


Menzell 


(1,1) 


3 


2 


2 


2 -> 5 -> 10 -s- 14 



Algorithm Performance 



5. Conclusion. It is a challenge problem to solve the polynomial systems with 
singular solutions. Various symbolic-numeric methods have been proposed for refining 
an approximate singular solution to high accuracy [3J 2J [51 |HJ [TU1 HH H3 [23 [33 131] ■ 
The breadth one case root refinement has been studied in [H [TU1 H3] • In this paper, 
we show how to apply strategies in [18] to reduce the size of matrices appeared in [51136] 
to obtain a more efficient algorithm for refining an approximately known multiple root 
for this special case. We have proved the quadratic convergency of the new algorithm 
when the approximate solution is close to the isolated exact singular solution. We also 
notice that when the singular solution x e is not well separated from other solutions 
of F, it is difficult to ensure that the approximate solution x will converge to x e . In 
[30] , they described an algorithm for computing verified error bounds for double roots 
of polynomial systems. We will explore ways of computing the certified bound for e 
to guarantee the convergency of our algorithm. It is also interesting to see whether 
the approach in the paper can be generalized to refine singular solutions when the 
Jacobian matrix is not of corank one. 
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Appendix. In the following, we suppose that there are no differential operators D(i, 0, . . . , 0) 
for i < k contained in Lu, otherwise, we can reduce them by Li. Here and hereafter, we 
always assume the coefficient of D(k, 0, . . . , 0) is one. 

Proposition 5.1. The formula \3.4 -?| ) is true for 1 < k < /i — 1. 

Proof. From (|3.35|l . we know that (|3.4ip is true for k = 1. Now let us assume (13.41 p 
is true for k. Let C^+j and C D k ^ denote the coefficients of D(a) in Pk+i and Qk+i 



respectively. In order to prove (|3.41[) for k + 1, we show that C^^ 1 = (k + l)cjj:+| by using 
the following relations repeatedly: 



C D(t) - (£>(<*))> Ql ^ °> 



(5.1) 



Let D(a) = D(ai, . . . ,a n ), denote ji = ■ ■ ■ = j ai = 1, • • ■ , j\ a \- ai +i = ■■■ = j\ a \ = i, 
where a, is the last nonzero entry in a, e.g., j\ — j% = jz = = 35 ~ 2,je = 3 for 

D(a) — D(3, 2, 1, 0). Since all derivatives in Qk+i and Pk+i are of order at least 2, we can 
start with \a\ = 2. 

1. If ai = and \a\ = 2, then we have 

C D(a) = 2a 2jl a fc-l,J2 H 1" ( k ~ l) a fc-l,Ji a 2,j 2 

+ 2o2, J - a Ofc_i lJ - 1 + ■ ■ ■ + (k - l)ak-ij 2 a 2 , J1 
=(k+ l)(a2,j 1 a k -i,j 2 + H0fe_i )J - 1 02,j a ) 

2. If ai and \a\ = 2, 

C^:, 1 = a kdl +ka k , 31 =(k + l)a Ml = (k + lJcQ+J. 

3. If ai ^ and \a\ > 2, by induction, 

(£>(«)) - /to * 1 (I3( Q )) 

n 

- ("i ~ l K\^ a)) + E °y (2« 2 „<% 2 ( D ( Q )) + - + (fc - l)a*-i,i^ #j(fl(a)) ) • 

j=2 
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While based on (|5.ip . we have C^+| = C^K D ,^ and 

n 

n 

= QlC, *|(D( Q) ) ( 2a2 J C '*xi 3 2 (I5(a)) +•'• + (*- 1 ) a fe-l J C* 1 1 * 3 (I5( Q ))) 



3=2 



= (* + i)Ci? W a)) = (* + i)cgcir 

4. If Qi = and |q| > 2, we have 

C D(a) -^Jl^ (£>(«)) + ■ ■ ■ + ^fc.JlC^fDfo 



H h 2a2, J |„ | C" <i) ^| (D(a)) H h fc • Cj 1 ^ (£>(<*»• 



By (E3), we have (7^ flW) = 02, J -C^* {i J {a)) + ■ ■ ■ + a*-ij^ (DW) . 
For 2 < p < fc, we collect all items with respect to Cp,^ : 

P%>,J1 (D(a)) + a PJ'i (^SjaCs'l^pta)) I" ~ P) a k-p,j 2 C 9 ^ 4- 2 (D(a))) 

H h a f>J'l ( 2a2 ' J l°l C '*j 1 *j|,!| (0(a)) + K fc ~-P) ' afe -P.J|c| (:7 * 3 1 1 * j|c>| (I3( a ))^) 

= P ■ a P,il (15(a)) + ( fc - P + l) a P,il (D(a)) = ( fc + !) a P,Jl C*^ (15(a))- 

Hence, = (fc + 1) Ep=2 a^C^i)) = ( fc + l)c£$ • 

Hence, (|3.4ip is true for 1 < fc < fi — 1. □ 

Proposition 5.2. In the proof of Theorem \3.1S\ we claim that F^(Li) is a linear- 
combination of {La , L\ , , Li } . 

Proof. For i = 2, since r^(Li) = -L\, we can reduce r^(L2) by Li to a differential 
operator which does not contain 73(1, 0, . . . , 0), denoted by f ^(£2). Since L2 is constructed 
by Theorem [2l] from L = 73(0, . . . , 0) and Li = 73(1, 0, . . . , 0), 73(2, . . . , 0) is the only 
second order derivative contained in L2 with coefficient one [181 Lemma 3.3]. Moreover, by 
Theorem l3.3l the coefficient of 75(2, . . . , 0) in T^(L2) is ^j. Since P^j(L2) is an approximate 
basis, due to the closedness, we have 

*i(fa(Z a )) = 

Therefore, we have 

Fh(7/ 2 ) = -rL 2 . 
a 2 

Hence, F^Z^) is a linear combination of {Lo, L\ , L2}. Let us assume that for all < i < fc— 1, 
the proposition is true. We can reduce F^(Lfc) by {Lo, Li , . . . , Lk-i} to a differential operator 
which does not contain differential operators D(i, 0, . . . , 0), denoted by F^(Lfc). Since F^(Lfc) 
is an approximate basis, due to the closedness, we have 

*i(r ft (I k )) = 4 £ *-i> 

a K 

therefore, r^(Lfc) — \Lk. Hence, T^(Lfe) is a linear combination of {Lo, Li, . . . , Lfc}. □ 



